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ABSTRACT 

Star formation plays an important role in the fate of interacting galaxies. To date, 
most galactic simulations including star formation have used a density-dependent star 
formation rule designed to approximate a Schmidt law. Here, I present a new star 
formation rule which is governed by the local rate of energy dissipation in shocks. The 
new and old rules are compared using self-consistent simulations of NGC4676; shock- 
induced star formation provides a better match to the observations of this system. 

Key words: galaxies: mergers - galaxies: kinematics and dynamics - galaxies: struc- 
ture 



1 INTRODUCTION 

Numerical simulations of galaxy formation and interactions 
often include rules for star formation. In many cases, the 
local rate of star formation, p*, is related to the local gas 
density, p g , by a power law: 

p, = C*pl, (1) 

where C» is a constant. This prescription has been justified 
both empirically and theoretically. From a n empirical per - 
spective, Q resembles the "Schmidt law" dSchmidd fl959'). 
which relates star formation per unit surface area, E* to gas 
surfac e density, S g ; a recent determination of the Schmidt 
law feennicutdli99a) is 

M / v \ 1 ' 4±0 ' 15 
E, = (2.5 ± 0.7) x 1(T 4 ( ,„ Sg , ] . (2) 

V ; yrkpc 2 \lM Q pc- 2 J y ' 

iMihos. Richstone fc Bothunl (ll99lT) adopted a rule equiva- 
lent to Q with n = 2 as an approximation to a Sch midt 
law with index ~ 1.8, while lMihos fc Hernauisd lll994aT ) took 
n = 1.5, and presented numerical tests showing this gave a 
reasonable match to a Schmidt la w with index ~ 1.5 . A more 
theore tical approach, adopted bv iKat A J1992D and lSprineell 
( 2000), sets p, = Pg/t*, where t*, the time-scale for star for- 
mation, is basically proportional to the local collapse time 
of the gas, (Gp g ) -1 ^ 2 . This yields Q with n = 1.5 and 
G» = G 1/2 C'«, where C, is a dimensionless constant. 

Despite the apparent convergence of observation and 
theory on the index n = 1.5, it's unlikely that Q really 
captures the process of star formation. For one thing, only 
about 1 per cent of the gas actually forms stars per collapse 
time (Gp g ) -1//2 ; in other words, consistency with the ob- 
servations implies that C* ~ 10 -2 . Current theories of star 
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formation don't offer any straightforward way to calculate 
this quantity. Gravitational collapse is evidently not the lim- 
iting factor which sets the rate of star formation; "feedback" 
processes are important in determining the fraction of avail- 
able interstellar material which ultimately becomes stars. 
While several groups have now devised simulations includ- 
ing s tar formation regulated by various fo rms of feedback 
(e.g. iGerritsen fc Ickdll997t ISpringej EoOOl. this approach 
still has some way to go. 

Moreover, models based on Q don't reproduce the 
larg e-scale star formation seen i n many interacting galax- 
ies. IMihos. Richstone fc Both un ( 1992) found that most 
of the star formation was con fined to the central regions 
of the ir model galaxies, while IMihos. Bothun fc Richstone! 
( 1993) noted that Q underestimated the rate of star for- 
mation in regions where interstellar gas exhibits large veloc- 
ity dispersions and gradients, or where the gas appears t o 
be undergoing strong sho cks. IMihos fc Hernq uisd (^994b), 
IMihos fc Hernauisd (jl996l> . and ISpringell J200(J) have mod- 
eled the central bursts of star formation seen in ultralu- 
minous infrared galaxies llSanders fc Mirabe 1 I1996L and ref- 
erences therein). But nuclear starbursts, while a necessary 
stage in the trans formation of merger remnan ts into ellip- 
tical galaxies fe.g. lKormendv fc Sanderall992t) . may not be 
sufficient to accomplish this transformation. For example, 
merger-in duced starbursts create mass ive young star clus- 
ters (e.g. Wh itmore fc Schweizerl fl995) which may subse- 
quently evolve into globular clusters, but these clusters will 
be confined to the nuclei of remnants unless the starbursts 
are spatially extended. 

In view of these considerations, it's worth examin- 
ing alternatives to Q. One long standing idea is that 
collisi ons of molecular clouds tr igger of star formation 
fe.g. lScoville. Sanders, fc Clemenslll986lh This has been im- 
plemented in "sticky particle" schemes which model the in- 
terstellar medium as a collection of discrete clouds under- 
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going inelastic collisions (e.g. lOlson fc Kwan|[l99(j iNoguchil 
1991). However, molecular clouds have relatively long mean 
free paths, and direct collisions between clouds at velocities 
of 10 2 km/s or more may resul t in d isruption rather than 
star formation. Ijog fc Solomon! dl992h proposed a model of 
shock-induced star formation in interacting galaxies; specif- 
ically, they suggested that fast collisions between extended 
HI clouds create a high-pressure medium which compresses 
pre-existing molecular clouds and thereby induces bursts of 
star formation. High-pressure regions, and especially large- 
scale shocks in colliding gala xies, may favor the fo rmation 
of bound star clusters l)Elmegreen fc Efremovlll997l l. 

Recent observations at optical, infrared, and radio 
wavelengths continue to reveal large-scale star forma- 
tion in interacting galax ies, including NGC 4038/9, 
Arp 299, and NGC 4676 iWhitniore fc Schweizel Il995t 
Vigroux et alJ Il996t iMirabel et all Il99& Iffibbard fc YurJ 
19991: lAlonso-Herrero et alJ 1200(1 Ixu et all l2000t 
de Griis et alJ I2003T) . Several of these studies explicitly 
invoke shock-induced star formation in discussing the 
observational data. In this paper I focus on NGC 4676, a 
strong ly interacting pair of dis k galaxies with lo ng tidal 
tails (iToomre fc Toomrd Il97l hereafter TT72; iToomrd 
Il977l) . as a test-case for models of shock-induced and 
density-dependent star formation. § 2 describes the star 
formation algorithms. § 3 presents simulations of NGC 4676, 
and contrasts the results of shock-induced and density- 
dependent star formation rules. Conclusions are given in 
§4. 



2 SIMULATING STAR FORMATION 

Within the fram ework of "Smoothed Particle Hydrody- 
namics" or SPH iLucvlll977t iGingold fc MonaeharJ fl977t 
iMonaghanl Il992l and references therein), the interstellar 
material can be modeled as an isothermal gas with a 
representative temperature T ~ 10 4 K. This approach is 
less com plex than others which i nclude radiative cool- 
ing (e .g. iBar nes fc Hernquistl Il996l) . heating by massive 
stars ([Gerritsen fc I ckc 1997), or other feedback effects 
( Springcl 2000). However, previous studies indicate that sim- 
ple isothermal models capture much of the relevant dynami- 
cal behavior of interstellar gas in interacting galaxies, includ- 
ing large-scale shocks, nuclear inflows, and the formation of 
extended gas disks in merger remnants ([B arnes fc Hernquistl 
Il996l: iMihos fc Hernauistlll996l : lBarnedl2002ll 

With the gas temperature fixed, there remain two inter- 
esting local quantities which are easily evaluated in an SPH 
simulation: the gas density, p g , and the rate of mechanical 
heating due to shocks and PdV work, ii. (An isothermal 
model assumes that radiative processes locally balance me- 
chanical heating or cooling.) In galactic simulations, the gas 
temperature is much lower than the virial temperature, so 
most shocks are highly supersonic and shock heating ex- 
ceeds PdV work by about two orders of magnitude. Thus 
shocks can be identified as regions with ii significantl y larger 
than zero; this criterion serves to locate shock fronts jBarnesI 
I2002H 

Making use of both p g and ii, I consider a general star- 
formation prescription of the form 

p* = C, pg MAX(u, 0) m , (3) 



where the MAX function, which returns the larger of its two 
arguments, is used to handle situations in which u < 0. This 
prescription is implemented in a probabilistic fashion; the 
chance of gas particle i undergoing star formation in a time 
interval At is 

Pi = apf l ~ 1) MAX(u i ,0) m At, (4) 

where pi and iii, the gas density and heating rate for par- 
ticle i, are defined in Appendix A. At each time-step, p,: is 
computed for each gas particle and compared with a ran- 
dom number X freshly drawn from a uniform distribution 
between and 1. If pi > X then the gas particle is instantly 
and completely converted to a star particle, which subse- 
quently evolves collisionlessly. For later analysis, each star 
particle formed during a simulation is tagged with the value 
of the time t at its moment of birth. 

While © can encompass a wide range of rules, this 
paper focuses on two limiting cases. First, setting n > 1 
and m — yields density-dependent star formation, in 
which pi oc p'" -1 ' with no dependence on u, and regions 
of higher density produce more stars. Second, setting n = 1 
and m > yields shock-induced star formation, in which 
Pi oc MAX(?i,0) m ; star formation is strictly limited to re- 
gions with u > 0, which typically correspond to shocks. 
More general rules with both n > 1 and m > are beyond 
the scope of this paper, but it may be possible to anticipate 
their behavior using experience gained with the simple cases 
described here. 

Simulations of shock-induced star formation are subject 
to resolution effects. The width of a shock region in SPH is 
~ 3 times the local smoothing scale hi iMonaghanl I1992T) 
(in these experiments, hi is set to enclose A?" sm ooth = 40 
gas particles within a radius 2hi of each gas particle i; see 
Appendix B). Consider a situation in which a gas particle 
encounters a shock of some fixed strength. The time that 
the particle spends in the shock scales as hi, while the to- 
tal amount of energy it dissipates is independent of hi, so 
iti oc hj 1 . Thus the net probability that the gas particle is 
transformed into a star particle while passing through the 
shock is proportional to h^~ m \ The specific choice m = 1 
should produce results which are approximately independent 
of the spatial resolution. For other choices for m the star 
formation rate will depend on resolution; to obtain roughly 
equivalent results in simulations with different spatial reso- 
lutions, it will be necessary to adjust the parameter C*. 

Resolution effects also arise in simulations with density- 
dependent star formation rules. Such simulations do a good 
job of reproducing the observed r elationship between gas 
density and st ar formation rate jMihos fc Hernauisdir994al : 
Springel 2000). But current SPH simulations generally don't 
resolve the vertical structure of gas in galactic disks, so 
this success may be due, in part, to compensating errors 
at different distances from the disk midplane. More criti- 
cally, although shocks with Mach numbers of M ~ 10 to 30 
will boost the de nsity of a n isothermal gas by factors of 
M 2 ~ 10 2 to 10 3 JShulll992l) . these high-density regions are 
not resolved in existing SPH calculations, and the localized 
bursts of star formation which should result have not been 
seen in the numerical models. In principle, simulations us- 
ing density-dependent rules can exhibit shock-induced star 
formation, but the computing power required to achieve the 
necessary resolution may not yet be available. 
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Figure 1. The Mice, NGC4676, and the model. North is up, west is right . The north-most hulk is NGC 4676a, while its companion to 
the southeast is NGC 4676b. Far left: deep optical image and HI contours iHibbard fc van Gorkomlll996l) . Left: optical image, stretched 
to show the bodies and tails. Right: self-consistent N-body model. Far right: declination versus line-of-sight velocity plane; contours are 
HI data, points are N-body model. 



3 MODELING THE MICE 

On account of their rounded bodies and long tails, the two 
galaxies making up NGC 4676 were du bbed the "Playing 
Mice" bv IVorontsov-Vervaminovl lll958l) . The far left panel 
of Fig. ^ shows a deep optical image , overlain by contours 
of HI (Hibbard & van Gorkoml H996l. while the adjacent 
panel shows a shallower image. The individual galaxies are 
NGC 4676a, to the north, and NGC 4676b, to the southeast. 
The length and straightness of the bright northern tail imply 
that we are viewing the aftermath of a nearly-direct passage 
roughly edge-on to the orbit and disk planes, while the curve 
of the fainter southern tail is consistent with a more face- 
on view of an inclined disk. The fairly equal lengths of the 
two tails suggests that the galaxies involved had roughly 
equal masses. Along with long-slit da ta indicating that both 
"hulks" rotate with north receding JpTurbidge fc Burbided 
ll96ll : lThevs. Spiegel, fc Toomrejll972f) . these considerations 
led TT72 to a simple model for this system. Their model re- 
quired that the northern hulk recede faster than the south- 
er n one, a predicti on soon confirmed with Ha spectroscopy 
bv lStocktonl lll974tl . Most kinematic studies, whether b ased 
on Ha spectros copy iM ihos 1 Bothun fc Richsto ndll993}) . HI 
interferometry (Hibbard & van Gorkoml ll99fE " or CO in- 
terferome trv jYun fc Hibbard! 1200 lT) . yield results consis- 
tent with Stoc ktorjs. On the other hand, Ha spectroscopy 
bv lSotnikova fc Reshetnikovl jl99ctl indicated a remarkably 
large difference of ~ 300 km s -1 between the nucleus of 
NGC 4676a and its associated tail, and a large velocity dis- 
persion within the tail its elf; it's tempting to interpret this 
as measurement error, but lThevs. Spiegel, fc Toomrd <ll972T) 
reported somewhat similar results, so the situation is not en- 
tirely clear. Nonetheless, HI interferometry appears the best 



choice for tracing the large-scale dynamics of the system, 
inasmuch as it provides a full map of the velocity field, and 
neutral hydrogen is unlikely to be disturbed by star forma- 
tion. 

3.1 Dynamical model 

The model of NGC 4676 used in this paper is based on work 
with J. Hib bard (Hibbard & Barnes, in preparation; see also 
iBarnes! 11998). In view of some concerns noted below, this 
model is somewhat preliminary, but it seems to describe the 
overall evolution of the system reasonably well. 

The initial conditions for the individual galaxies each 
have three com ponents: a bulge with a shallow cusp 
jHernauistlll990|). an exponential d isk with constant scale 
height jFreemanll970l;ISpitzeJl'942T) . and a dark halo with a 
constant-density core, base d on a "gamma model" jDehnenl 
119931: iTremaine et alJll994h . Density profiles for these com- 
ponents are 

Pbuigc oc r _1 (r + a b ui gc )~ 3 , 

Pdisk oc exp(-i?/i?disk) sech 2 (2:/zdisk) , (5) 
Phalo cc (r + a ha i )~ 4 , 

where r is the spherical radius, R is the cylindrical radius 
in the disk plane, and z is the distance from the disk plane. 
Numerical results are presented using a system of units with 
G — 1. In these units, the components have masses Mb u i gc = 
1/16, Af disk = 3/16, M ha io = 1, and length scales a bu i g c = 
0.04168, R disk = 1/12, 2 disk = 0.005, and a ha i G = 0.1. With 
these parameters, the disk has a rotation period of ~ 1.0 
time units at R = 3J?disk- Scale factors relating simulation 
units to real physical quantities will be determined as part 
of the model-matching process. 
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Beginning with TT72's model, we ran experiments sys- 
tematically varying the relative orbit and orientation of the 
two identical galaxies. The following parameters produced 
the reasonably good fit to NGC4676 shown on the right 
and far right of Fig.0 First, the initial orbit was parabolic, 
with a pericentric separation r p = 0.25 = 3i?disk- Second, 
the disks of NGC 4676a and NGC 4676b were initially in- 
clined to the orbital plane by i n — 25° and ib = 40°. Third, 
the initial angle between each disk's line of nodes and the 
direction of separation at pericenter was cj a = —30° and 
cob = 60°. Fourth, the system was evolved 1 time unit past 
pericenter. We used an interactive display program to ro- 
tate the model to best match the observed morphology and 
kinematics NGC 4676; further details of the matching pro- 
cess will appear in a subsequent paper. 

Comparing the model with the observations, the main 
discrepancy is the position angle of NGC 4676b's bar; optical 
images show the bar runs roughly NNE-SSW, while the bar 
in the model is roughly N-S or even NNW-SSE. The pat- 
tern speed of a bar depends on the details of the rotation 
curve, so this discrepancy might be resolved by adjusting 
the mass model used for NGC 4676b. It's also possible that 
NGC4676b's bar predates its encounter with NGC 4676a, 
in which case the bar's cur rent position angle would depend 
on it s phase at pericenter llGerin. Combes, fc AthanassoulJ 
1990). However, the hypothesis that NGC4676b's bar was 
induced by the tidal interaction is more parsimonious and 
has not yet been tested by experiments using a range of mass 
models. 

The time yielding the best match between the model 
and NGC 4676 proved somewhat ambiguous; by adjust- 
ing the viewing angle and scaling factors, acceptable 
matches to the tail morphology and kinematics were ob- 
tained for t ~ 0.75 to 1.0 time units after pericenter. 
At t = 0.875, a good match was obtained by equat- 
ing one length unit to 80" = 35.7 kpc, and one veloc- 
ity unit to 180 kms" 1 = 0.184kpcMyr _1 . This scaling as- 
sumes that NGC 4676's distance is 92 Mpc, consistent with 
H = 72kms~ 1 Mpc" 1 and NGC 4676's syste mic velocity 
of 6600 km s" 1 jHibbard fe van Gorkoml 119961) . One time 
unit then works out to 194 Myr, and pericenter occurred 
~ 170 Myr ago. Finally, one mass unit is 2.69 x 10 11 M . 

Given this scaling, the initial parameters of the model 
galaxies can be computed. For example, the initial disks 
have reasonable scale lengths -Rdisk — 3.0 kpc and masses 
Mdisk — 5.0 x 10 10 Mq; on the other hand, they have maxi- 
mum circular velocities of ~ 290kms _1 , which seem rather 
high and might be repaired by increasing the halo scale ra- 
dius Ohalo- 

Since TT72's work, NGC 4676 has been modeled by 
several other groups. iMihos. Bothun fc Richstoriel lll99St) es- 
sentially produced a self-consistent clone of TT72's model, 
even adopting a short-period elliptical orbit like the one 
TT72's test-particle methodology dictated; they favored 
an earlier viewing time which matched the position an- 
gle of NGC 4676b 's bar b ut produced tails s ignificantly 
shorter th a n thos e observed. iGilbert fc Scllwood ( 1993|) used 
IStocktonl 's lll974h velocity data to construct a self-consistent 
model with a parabolic initial orbit. In most respects, their 
model is very similar to the one described here, although we 
favor a closer passage to i ncreas e the heft of the tidal tails. 
ISotnikova fc Reshetnikovl jl998T) produced a test-particle 
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Figure 2. Global star formation in simulations of NGC 4676. The 
top panel shows M* , the rate of star formation, while the bottom 
panel shows the mass of gas remaining; both quantities are given 
in simulation units (G = 1). Pericenter occurred at t = 0; the hor- 
izontal bar labeled "now" indicates the range of times matching 
the morphology and kinematics of NGC 4676. Solid line: density- 
dependent star formation with n = 1.5, m = 0, C* = 0.025. 
Broken lines: shock-induced star formation; dashed is n = 1, 
m = 0.5, C, = 0.5, dotted is n = 1, m = 1, C* = 0.25. Note 
that the shock-induced models produce strong bursts of star for- 
mation at pericenter (t = 0), while the density-dependent model 
exhibits little or no immediate response. 



model, restricted to NGC 4676a, which matched the general 
trend of their velocity data but not the large velocity dis- 
persion they reported within the northern tail. Their model 
requires relatively massive dark halos to account for these 
high tail velocities. 

3.2 Interaction-induced star formation 

Fig- El compares global star formation rates for three differ- 
ent simulations of NGC 4676. To illustrate the range of pos- 
sibilities, I contrast one simulation with density-dependent 
star formation (solid line: n = 1.5, m = 0, C* = 0.025) and 
two simulations with shock-induced star formation (dashed 
line: n = 1, m = 0.5, C* = 0.5; dotted line: n = 1, m = 1, 
C» = 0.25). All of these simulations were based on the initial 
orbit and disk orientations described in § 3.1; the gas was 
initially distributed like the disk stars, and the gas mass in 
each disk was M gas = Mdi s k/8 = 3/128. Further numerical 
details are given in Appendix B. 

In simulations with density-dependent star formation, 
setting the constant C* = 0.025 yielded a baseline star for- 
mation rate per disk of M* ~ 0.00136 (see Appendix C); at 
this rate, an undisturbed disk would consume its gas on a 
time-scale of M gas /M* ~ 17 time units (~ 3.3 Gyr). In the 
simulation with shock-induced star formation and m — 0.5, 
a similar baseline rate of M* ~ 0.00126 was obtained by 
setting C* = 0.5; this gave reasonable behavior throughout 
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Figure 3. Dynamical evolution and star formation in simulations of NGC4676. Top row: old stellar distribution (contours) and gas 
density (halftone) at equally-spaced times between t = (pericenter) and 0.5. Middle row: old stellar distribution (halftone) and star 
formation (points) for the density-dependent model (ra = 1.5, m = 0, C* = 0.025). Bottom row: old stellar distribution (halftone) and 
star formation (points) for the shock-induced model (n = 1, m = 0.5, C* = 0.5). Red points have ages r < 1/16 time units, green have 
r < 1/4, and blue have r < 1. See http://www.ifa.hawaii.edu/~barnes/research/intcraction_modcls/mice/ for animated versions of this 
and other figures. 



the encounter. On the other hand, it was necessary to adopt 
a much lower baseline rate of M, ~ 0.000175 in the simu- 
lation of shock-induced star formation with m = 1 to avoid 
using up most of the gas in the earliest stages of the en- 
counter; setting C» = 0.25 insured that this case ultimately 
consumed about as much gas as the other two. 

As Fig. |2]shows, the density-dependent simulation pre- 
dict that the star formation rate remains fairly constant un- 
til some time after the first passage. Eventually, as the gas 
encountered shocks and gravitationally interacted with the 
tidally induced bars in the stellar disks, it began accumu- 
lating in the centers of the galaxies, driving a fairly gradual 
increase in star formation. The two shock-induced simula- 
tions, on the other hand, produced sharp bursts of star for- 
mation as the disks interpenetrated (t ~ 0), followed by brief 
lulls as the disks separated. Soon thereafter, star formation 
rates rose again as gas flows intersected in the perturbed 
disks; this produced the broad peaks, centered oni ~ 0.6, 
accounting for most of the stars formed before t = 2. As one 



might expect, the simulation with m = 1 exhibited the most 
extreme form of this behavior, with a higher peak at t = 
and a deeper minimum immediately thereafter. 

These various star formation histories have interest- 
ing observational implications. According to the density- 
dependent simulation, star formation in NGC 4676 has just 
begun responding to the encounter, and the peak star for- 
mation rate will occur several hundred Myr from now. In 
contrast, the shock-induced models predict that we are ob- 
serving the system a short time after the broad peaks of 
star formation at t ~ 0.6. This difference in star formation 
history is fairly insensitive to the choice of C», which influ- 
ences the overall amplitude of star formation but does little 
to change the timing of starbursts. 

Fig. [3] contrasts the spatial distribution of star forma- 
tion in density-dependent and shock-induced simulations. 
The top row of images set the stage by showing distribu- 
tions of old disk stars and gas, taken from a simulation 
without star formation. The middle row shows distributions 
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Figure 4. Simulations matching the present configuration of NGC 4676. Left: old stellar distribution (contours) and gas density (halftone). 
Middle: old stellar distribution (halftone) and star formation (points) for the density-dependent model (n = 1.5, m = 0, C* = 0.025). 
Right: old stellar distribution (halftone) and star formation (points) for the shock-induced model (n = 1, m = 0.5, C* = 0.5). Colors 
indicate ages as in Fig. l3l 



of old disk stars and star particles formed according to the 
density-dependent rule (n — 1.5, m = 0, C» = 0.025); colors 
indicate ages of star particles, with red being youngest and 
blue oldest. The bottom row presents results for a simula- 
tion with shock-induced star formation (n = 1, m = 0.5, 
C, = 0.5) using the same color coding. Results for the other 
shock-induced simulation are not shown since they are very 
similar. 

Comparison of the middle and bottom rows of Fig. [3] 
shows that the star formation rule clearly influences the dis- 
tribution of newly-formed stars. In the density-dependent 
model, star formation was always strongly concentrated to- 
ward the central regions of the disks. Some star formation 
did occur outside the centers, but it was relatively weak 
and short-lived. In contrast, the shock-induced models pro- 
duced large, spatially-extended bursts of star formation as 
the disks collided at t = 0, and relatively high levels of ac- 
tivity persisted in the bridge and tail regions until t ~ 0.4. 

By the last time (t = 0.5) shown in Fig. |3] products 
of star formation are much more widely distributed in the 
shock-induced simulations than in the density-dependent 
simulation. Rather interesting in this regard are the gas-poor 
features which extend to roughly five o'clock from the upper 
galaxy and to eleven o'clock from the lower galaxy. These 
regions were swept free of gas during the interpenetrating 
encounter of the two disks. In the simulations with shock- 
induced star formation, some of this gas was converted to 
stars before losing much momentum; as a result, star formed 
in the pericentric bursts are scattered throughout these gas- 
poor regions. By contrast, these regions are bereft of both 



gas and star-formation products in the density-dependent 
simulation. 

Fig. 2] presents views which approximate the observed 
morphology and kinematics of NGC 4676. These views are 
plotted at time t = 0.875, but fairly similar results are 
found between t — 0.75 and t = 1. The simulations 
shown here are the same ones shown in the last figure, 
with the density-dependent case in the middle and the 
shock-induced one on the right. In both cases, most of 
the ongoing star formation is now confined to central re- 
gions of the galaxies. However, products of star forma- 
tion are found throughout the tails and other tidal fea- 
tures, especially in the simulation with shock-induced star 
formation. These results invite comparis on with observa- 
tions of past star formation in NGC 4676. IStocktonl dl974l) 
found that the spectrum of NGC 4676 is dominated by A 
stars, and attributed this to "rapid, widespread star for- 
mation that effectively ceased at least 5 x 10 7 years ago"; 
this diagnosis is supported by more recent observations 
( Mihos 1 Bothun fc Richstonej|l993|: [Hibbard fc van Gorkoml 
Il996t ISotnikova fc Reshetnikovlll99£l Ide Griis et alJl2003l) . 
The shock-induced simulations produced just such bursts at 
first passage (t ~ 0), and the products of these bursts are 
widely distributed throughout the tidal tails at later times. 
In contrast, the density-dependent simulation had far less 
star formation outside the galactic centers, and the first pas- 
sage did not induce a prompt burst of extended star forma- 
tion. 

NGC 4676 has regions of ongoing sta r formation which 
are evident in Ha images jHibbardlll99l Fig. AIII-1). The 
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Table 1. Mass of stars formed since start of interaction, expressed 
as a percentage of initial gas mass. 



n 


m 


extended 


central 


total 


1.5 


0.0 


1.4 


9.2 


10.6 


1.0 


0.5 


5.9 


14.7 


20.6 


1.0 


1.0 


5.2 


13.6 


18.8 



central regions of the two galaxies a re quite prominent in Ha 
jMihos. Bothun fc Richstonelll993l) : central activity is seen 
in all the simulations. Outside of the centers of the galax- 
ies, star-forming regions are observed in the straight tail ex- 
tending to the north of NGC 4676a. None of the simulations 
presented here reproduce this phenomenon, although the 
shock-induced models come closer inasmuch as they show 
enhanced star formation in their tails until t ~ 0.4; the 
density-dependent model seems less promising in this regard 
since star formation in its tails declines after first passage. 
Ha images also show ongoing star formation at the base of 
the curved tail to northeast of NGC 4676b. This region has a 
relatively high surface bri ghtness and is consp icuously blue 
in the ACS images (e.g. Ide Griis et aljl2005tl . It is nicely 
reproduced in the shock-induced simulations; gas from the 
curved tail falls back into its original disk, and the resulting 
shock sustains a small, coherent patch of extra-nuclear star 
formation. 

On the whole, the simulations with shock-induced rules 
qualitatively match most of the recent large-scale star forma- 
tion in NGC 4676. Eventually, this comparison can be made 
more precise by modeling the spectrophotometric evolution 
of star particles, but such refinements are beyond the scope 
of this paper. The balance between central and extended star 
formation can be quantified by measuring the mass in stars 
formed since the start of the interaction. Specifically, the 
production of stars was integrated from t = —0.25, which is 
when Fig. [5] shows the first increase in interaction-induced 
star formation, until t = 0.875; particles then found at radii 
r > 0.15 from the centers of both galaxies were considered to 
be part of the extended tidal features. Table^list the results 
for all three experiments. At the epoch matching the present 
configuration of NGC 4676, the shock- induced models have 
formed about twice as many stars as the density-dependent 
model; of course, this ratio depends on the values chosen 
for C». In the shock-induced simulations, 28 percent of the 
newly-formed stars are found within the tidal features, while 
the corresponding figure for the density-dependent simula- 
tion is only 13 percent. These percentages should be insensi- 
tive to reasonable variations of C*, and preliminary experi- 
ments with other parameter values support this expectation. 

3.3 Starbursts in mergers 

The future evolution of the NGC 4676 system may seem a 
matter of conjecture. Nonetheless, it's worth asking when 
these galaxies are likely to merge and if their merger will 
produce an ultraluminous IR galaxy; even if we can't check 
the specific predictions for this system, the results may pro- 
vide insight into other merging galaxies. I therefore present 
some predictions for the fate of NGC 4676. 

Fig. extends the comparison of global star formation 
rates in NGC 4676 some 600 Myr into the future. In these 
simulations, star formation is modulated by two key events. 
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Figure 5. Global star formation in simulations of NGC 4676. 
First pericenter occurs at t = 0, second pericenter occurs at t ~ 
2.2, and the galaxies merge at t ~ 2.8. Line types are identical to 
those in Fig. [2] 

First, as a result of orbital decay, the galaxies have a second 
and much closer passage at t ~ 2.2; this passage has little 
immediate effect in the density-dependent model, but pro- 
duces definite bursts in the shock-induced models, especially 
in the one with the steeper m = 1 dependence on u. Second, 
the galaxies merge at t ~ 2.8, producing dramatic bursts of 
star formation in the density- dependent model and in the 
shock- induced model with m = 1, but only a modest burst 
in the shock-induced model with m — 0.5. 

In the shock-induced simulations, the bursts at the 
second passage were produced by the same mechanism 
responsible for the comparable, albeit somewhat larger, 
bursts at first pericenter: interpenetrating encounters be- 
tween two gassy disks. While much of the gas originally 
populating these disks was driven into the nuclei or con- 
verted into stars, gas retur ning from th e tidal tails set- 
tled back into the disks (e.g. lBarnesll2002T) ; thus, even after 
the fireworks attending their first passage, the disks had 
enough extended gas to produce brief starbursts as they 
came together for a second time. It's worth noting that 
such "second-passage" activity is apparently occurring in 
Arp299, whic h exhibits both an extende d starburst only 
4 Myr old (e.g. lAlonso-Herrero et alj|2 000[) and a lon g tidal 
tail with an age of ~ 750 Mvr (iHibbard fc Yunlll999l) . 

After the second passage, the orbits of the galaxies de- 
cayed very rapidly, and a final, almost head-on encounter 
resulted in merger at about t — 2.8. As noted above, this 
event produced fairly intense bursts of star formation in two 
of the three cases studied here. In the density-dependent 
simulation, this burst resulted from a sudden increase in 
gas density as the nuclei finally coalesced. Similar results 
have been seen in other numerical studies using d ensity- 
dependent star formation; iMihos fc Herncmisd il996f) report 
bursts with peak star formation rates as much as 10 2 times 
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Table 2. Log median radii of burst populations. 
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Table 3. Mass of burst populations, expressed as a percentage of 
initial gas mass. 
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the baseline. The present calculation yielded a more mod- 
est burst, peaking a factor of ~ 8 above the baseline; the 
lower peak in the present simulation may be due to con- 
sumption of gas at times t 1 to 2 while the galaxies were 
still separate. ISpringel ((2000) , combining density-dependent 
star formation with a sophisticated feedback scheme, ob- 
tained final bu rsts comparable to the one shown here. On 
the other hand. lMihos. Bothun fc Richs tonc ( 1993) reported 
only a "weak merger-induced" starburst in their model of 
NGC4676; while they blame this fizzle on the long delay 
between first passage and merger in their dynamical model, 
it's possible that the steeper n = 2 exponent of their star 
formation prescription had some role in consuming almost 
all available gas before the final encounter. 

Up until this point, the two simulations with shock- 
induced star formation have produced fairly similar results. 
Nonetheless, as Fig. 0] shows, their final bursts are quite 
different; the m = 0.5 model yields a small burst with a 
maximum M* ~ 0.0085, while the m = 1 model yields a 
sharply-peaked burst with a maximum M, ~ 0.0285. In 
principle, this difference could originate at earlier times - 
for example, the m — 0.5 model might yield a smaller final 
burst because it had already used up its gas - but a de- 
tailed examination of the simulations does not support this 
idea. At the start of their final encounters, both simulations 
have comparable amounts of high-density nuclear gas, and in 
both simulations these central regions collided almost head- 
on, producing comparable amounts of dissipation. Thus, the 
difference in burst strengths reflects the different behavior 
of the m = 0.5 and m = 1 rules under nearly identical 
circumstances: in both CEIS6S, £L moderate amount of gas en- 
countered a strong shock as the nuclei collided. The steeper 
dependence on u in the m = 1 case converted a large frac- 
tion of the shocked gas into stars, while in the m = 0.5 case 
the probability of star formation is less strongly enhanced 
and less gas was consumed. W hile the final bursts are no t 
as dramatic as those reported bv lMihos fc Hernauisd il996fl . 
it's possible that further simulations of shock-induced star 
formation with m = 1 could produce bursts large enough to 
power luminous IR galaxies. 

Tables Hand 01 summarize properties of the starbursts 
produced in these simulations. Here a burst is defined as 
all stars formed in a certain time interval; these intervals 
bracket various peaks of M* in Fig.^l For the shock- induced 
simulations, burst "Bi" corresponds to —0.20 ^ t ^ 0.12, 
"B 2 " to 0.28 < t ^ 0.93, "B 3 " to 2.10 < t < 2.34, and "B 4 " 
to 2.69 ^ t ^ 3.00. For the density-dependent simulation, 
"B 2 " corresponds to 0.80 t sC 2.34, and "B 4 " to 2.69 



t ^ 3.15. In effect, bursts Bi and B3 occur when the galaxies 
interpenetrate at first and second passages, respectively, B 2 
is the broad peak of star formation in the individual galaxies 
after their first passages, and B 4 is the final merger-induced 
starburst. The last column in each table, labeled "all", gives 
figures for all stars formed during these simulations (— 1 ^ 
is: 4). 

The radial distribution of burst populations in the fi- 
nal (t = 4) merger remnants is characterized in Table |21 
Here, the center of each remnant r cen t was defined by lo- 
cating the minimum of the potential well; next, the median 
of \ri — r con t| for all particles i in each burst population 
was computed; common logarithms of these median radii 
are listed. On the whole, the burst populations in the shock- 
induced simulations are ~ 3 times more extended than those 
in the density-dependent simulation; this ratio applies not 
only to the stars produced by bursts B 2 and B 4 but also 
to the entire population of stars formed during the simula- 
tions. The extended distribution of star formation products 
in the shock-induced models may provide a bet ter match to 
galaxies like NGC7252, where ISchweizerl ( 1990) found spec- 
tra dominated by A stars both within the nucleus and at a 
projected distance of 7kpc. In addition, the final B 4 popula- 
tions in the shock-induced simulations, while still compact, 
may be extended enough to moderate the cen tral "spikes" of 
young stars that lMihos fc Hernauisd Jl994ch found in their 
simulated merger remnants. 

Table 01 lists masses of stars formed. These masses are 
given as percentages of the total initial gas mass, 2M gas = 
0.046875. The two shock-induced simulations, by design, 
formed comparable amounts of stars, while the density- 
dependent simulation ultimately formed somewhat more. 
Burst B 2 , triggered within the individual galaxies after their 
first passage, yielded 41 percent of all stars formed in the 
density-dependent model, but only 24 percent in the shock- 
induced models. The final burst B 4 accounts for between 
8.2 percent and 23 percent of the total star formation. 

Earlier merger simulations including shock-induced star 
formation have largely been based on "discrete cloud" mod- 
els of the ISM, m aking a d i rect co mparison with the present 
results difficult. iNoeuchil dl99lT) presented cloud-collision 
rates for several mergers; assuming that each cloud colli- 
sion resulted in star formation, he predicted that merging 
galaxies undergo "repetitive starbursts" as the relative orbit 
of their nuclei decays. The present calculations also produce 
bursts of star formation associated with pericentric passages, 
but orbit decay in these f ully self-co nsistent calculations is 
much more abrupt than in iNoguchl 's experiments, possibly 
because the self-gravity of the gas is included. Thus, while 
there is a good deal of temporal structur e in the st ar forma- 
tion rate, the nearly-periodic starbursts iNoeuchil reported, 
which seem to reflect the gradual decay of elo ngated orbits 
in a n early-harmonic core, are not seen here. lOlson fc Kwanl 
( 1990) modeled the evolution of a spectrum of cloud masses 
during a merger, and reported a dramatic increase in the rate 
of disruptive cloud-cloud collisions within the inner 2 kpc of 
their model galaxies. Assuming that star formation is as- 
sociated with disruptive and glancing cloud-cloud collisions, 
they obtained peak luminosities and luminosity-to-gas ratios 
typical of luminous IR galaxies. 
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4 CONCLUSIONS 

Density-dependent and shock-induced models of star forma- 
tion yield qualitatively and quantitatively different results. 
These differences arise despite the fact that the gas density 
p g and the dissipation rate it are not completely indepen- 
dent variables; indeed, they can be explained in terms of the 
global relationship between p g and u. Dissipation is a nec- 
essary precursor for any significant increase in central gas 
density, since gas can only accumulate in galactic centers as 
a result of an irreversible process. The rate of increase in gas 
density closely tracks the ne t energy radiated away in shocks 
l|Barnes fc HernauisdllQQfiL Fig. 6); globally, u and dp e /dt 
are strongly correlated. Thus, models of shock-induced star 
formation can respond promptly to external disturbances, 
while little activity occurs in density-dependent models un- 
til p g has had time to build up. The larger spatial extent of 
star formation in shock-induced models is, in part, a corol- 
lary of the earlier onset of activity in such models, since the 
gas is more widely distributed at earlier times. 

Density-dependent rules may be used to implement 
unified models of star formation in normal and interact- 
ing galaxies. The disks of normal galaxies and the cen- 
tral r egions of s t arbur st galaxies fit onto the same power 
law ijKennicuttl Il99ftf) . Thus, by setting C, to match 
the baseline rate o f star formation in unperturbed disks, 
iMihos fc Hernauistl (|l99fj) were able to produce starbursts 
compara ble to those inferred in ultraluminous infrared 
galaxies dSanders fc Mirabell Il996l) ; this represents a real 
success for the model defined by £Q| with n = f .5. 

But there is abundant evidence for large-scale star 
formation in interacting galaxies. This includes ongoing 
star formation in the "overlap regions" of systems like 
NGC 4038/9 and Arp299, as well as the Ha emission from 
the tails of NGC 4676. It also includes the A-star spectra 
seen in the tails of NGC 4676 and throughout the body 
of NGC 7252. Density-dependent star formation can't eas- 
ily explain these observations; interactions funnel gas from 
the disks to the central regions of galaxies, thereby promot- 
ing rapid star formation in galactic nuclei but reducing the 
supply of gas needed for star formation elsewhere. More- 
over, violent relaxation ceases l ong before b i nding energie s 
are effectively randomized (e.g. IWhitelll979t lBarneslll992|) . 
so merging is ineffective at transporting the products of star 
formation outward from nuclei to the bodies of merger rem- 
nants. Thus, density- dependent rules offer at best an incom- 
plete description of star formation in interacting galaxies. 

In simulations with shock-induced rules, the star for- 
mation rate depends on dynamical circumstances. Shocks in 
unperturbed disks are associated with transient spiral pat- 
terns, and matching the baseline level of star formation de- 
pends on reproducing the "right" level of spiral structure. 
This, in itself, is a tough problem, and it's not clear that 
shock-induced models will soon yield a unified description 
of star formation in normal and interacting galaxies. Until 
this becomes possible, unperturbed disks can't be used to 
set C», so the amplitude of a shock-induced starburst can't 
be predicted a priori. Nonetheless, shock-induced star for- 
mation is an important element of galactic collisions, and 
its implementation in numerical simulations is a useful step 
toward increased realism. 

These simulations have contrasted two limits of ©: 



density-dependent star formation, with n > I and m = 0, 
and shock-induced star formation, with n — 1 and m > 0. 
These cases were chosen as instructive examples, but 
is general enough to accommodate additional possibilities. 
First, setting n > 1 and m just slightly larger than would 
yield a modified law in which star formation is proportional 
to Pg but occurs only in regio n s wit h ii > 0; this resem- 
bles the rule adopted by iKatzl (I1992T) , who basically took 
p* oc py 2 but restricted star formation to regions with con- 
vergent flows. Notice that cases with < m < 1 are not 
continuous with the case m = 0; in the latter, star formation 
is strictly independent of u. Second, setting n > 1 and m > 
would yield hybrid rules in which star formation depends on 
both p g and u. The consequences of such rules can some- 
times be inferred from the limiting cases considered here. 
For example, in the simulations of NGC 4676 with shock- 
induced star formation, the gas involved in burst B3 has 
densities ~ 10 2 times higher than the gas involved in burst 
Bi; setting n > 1 would boost B3 relative to Bi by a factor 
of ~ lO 2 *-™ -1 - 1 , assuming that neither burst was limited by 
the supply of gas. 

Perhaps most exciting are the new avenues for research 
created by an alternative description of star formation in 
galaxy interactions: 

• In contrast to density-dependent star formation, which 
seem s fairly insensitive to mos t details of galactic encoun- 
ters jMihos fc Hernauistlll996l . different encounter geome- 
tries may yield a variety of shock-induced star formation 
histories. A survey of different encounters, along the line s 
of other surveys without star formation fe.g. lBarnesl [2002 % ). 
might establish if widespread star formation at first passage 
is generic or limited to a subset of close encounters, and also 
determine if shock-induced star formation can account for 
ultraluminous infrared galaxies. 

• Shock-induced star formation makes definite predic- 
tions about the timing of starbursts triggered by galactic en- 
counters. As already mentioned, models including the photo- 
metric evolution of starburst populations could sharpen the 
comparison between simulations and observations of systems 
like NGC 4676. In addition, these models might predict the 
ages of embedded young star clusters; such predictions could 
be checked by multi-object spectroscopy of the clusters. 

• In the context of merger simulations, shock-induced star 
formation also yields predictions for the spatial distribution 
and kinematics of starburst populations. It would be quite 
interesting to compare these predictions with observations of 
the distribution and kinematics of metal-rich globular clus- 
ters in elliptical galaxies, whic h may have been form ed in 
merger-induced starbursts fe.g. IZepf fc Ashmanlll993h . 

• Encounters between gas-rich galaxies at redshifts of a 
few seem a natural setting for shock-induced star forma- 
tion. Whi le many high-z galaxies h ave peculiar morpholo- 
gies (e.g. Ivan den Bergh et al.lll99rl) . the bridges and tails 
characteristic of low-z encounters (TT72) are not very evi- 
dent. Instead, the optical morphology of these objects may 
be dominated by rest-frame UV from widespread starbursts; 
if so, the photometric modeling approach described above 
could help interpret existing and forthcoming observations 
of interacting galaxies at intermediate and high redshifts. 
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APPENDIX A: SPH EQUATIONS 

I use an adaptive SPH code in which each gas particle i is 
assigned smoothing length hi which depends on local con- 
ditions. The gas density at p article i is estimated usin g the 
symmetric form suggested bv lHernauist fc Katzl |l989), 



(Al) 



where nrij is the particle mass and r\j = |r< — Here, 
Wij{r) = \[W{r,hi) + W(r,hj)], and W(r,h) is th e spline 
interpolation kernel iMonaehan fc Lattanziol ll985T> . which 
vanishes identically for r > 2h. 

As discussed in § 2, the code uses an isothermal equation 
of state, so the internal energy u, of each gas particle is con- 
stant. One may nonetheless compute the mechanical heating 
rate in associated with the pressure and (artificial) viscous 
forces on each particle i, and assume that this heating is 
exactly balanced by radiative processes. This requires only 
a modest number of computations over and above those re- 
quired to compute the forces themselves, and allows a global 
check of energy conservation. The mechanical heating rate 
is estimated using 



V; - V 



_d_ 

dri 
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where Pi = c\pi is the pressure at particle i for an isothermal 
gas with sound speed c s . Here, Hij approximates a combina- 
tion of bulk and von Neuman-Richtmyer artificial viscosity 
<Monaghadll992ft . 



n y = i (p* + p.)/2 
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otherwise 



where a and (3 are parameters of order unity, 

(r; -Tj) ■ (Vj -Vj) 
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hij = \{hi + hj), and r\ is a parameter of order 10 



APPENDIX B: SIMULATION DETAILS 

The disk galaxies used in the star formation simulations 
had a total of 87040 particles each, with 24576 for the gas, 
21504 for the stellar disk, 8192 for the bulge, and 32768 
for the halo. The gas was initially distributed like the stel- 
lar disk, and accounted for 12.5 per cent of the total disk 
mass. In simulation units, the gas had a fixed sound speed 
c s = 0.0966, mimicking a warm ISM with T ~ 10 4 K. 

These calculations were run with a SPH code featuring 
adaptive smoothing and individual particle timesteps. The 
smoothing radius hi of each gas particle was set so that a 



sphere of radius 2hi contains exactly N suloo th = 40 gas parti- 
cles. A modified Courant condition with C = 0.25 was used 
to determine individual particle timesteps, which in some 
cases were as small as At; = 1/32768 time units. The arti- 
ficial viscosity parameters appearing in 1A3I were a — 1.0 
and /3 = 2.0. Gravitational forces were calc ulated using a 
modified tree code with Plummer smoothing iAarsetbll963h 
and a smoothing length of e = 0.0125. Collisionless parti- 
cles were integrated with a fixed time-step of At = 1/512 
time units. Including the energy radiated (or absorbed) by 
the gas in order to maintain a constant temperature, energy 
was conserved to about 0.5 percent through the final merger 
of the galaxy models. Depending on the star formation pa- 
rameters, each simulation took between 126 and 245 hr on 
a 1 GHz processor; the run without star formation was the 
slowest since its gas attained the highest densities. 

In practice, the limited resolution of computer simula- 
tions introduces effects which complicate the interpretation 
of numerical experiments. Smoothing over a sphere of ra- 
dius 2hi inevitably suppresses structures on smaller scales. 
The gas disks used in this paper provide good examples; 
for the simulation parameters given above, only half of the 
disk particles have smoothing radii hi less than the verti- 
cal scale height ZdiskJ Clearly, the vertical structure of these 
disks is not resolved - and the same is true for most SPH 
simulations of interacting galaxies published to date. In this 
circumstance, SPH underestimates p g close to the midplane, 
and overestimates p s far from the midplane. 



APPENDIX C: ISOLATED DISK MODELS 

To better characterize the behavior of the galaxy models 
used in the simulations of NGC 4676, I ran some test cal- 
culations of isolated disks. These can be used to establish 
the baseline rate of star formation in non-interacting disks, 
and to examine the radial distribution of star formation for 
various choices of n and m. 

Fig. IC1I shows global star formation rates for isolated 
disk models like the ones used in the NGC 4676 simula- 
tions. As in Figs. [5] and the solid lines show results for 
the density-dependent rule (n — 1.5, m — 0, C* = 0.025) 
while the broken lines show results for shock-induced rules 
(dashed line: n = 1, m = 0.5, C* = 0.5; dotted line: n — 1, 
m = 1, C* = 0.25). Since the last of these yields a very 
low baseline rate, an additional calculation was run with C» 
increased by a factor of 8 (dot-dashed line: n = 1, m = 1, 
C* = 2.0). In every case, the isolated disk models proved 
stable for the duration of the calculation. Apart from the 
brief initial bursts in the shock-induced models due to im- 
perfect initialization of gas velocities, the only trends seen 
are slight declines in global star formation due to gradual 
depletion of the gas. 

Given the rather stable behavior of the isolated disk 
models, it's reasonable to measure radial star formation pro- 
files by counting stars formed in concentric annuli. Fig. IC2I 
plots E* as a function of cylindrical radius. Much as ex- 
pected, the density-dependent rule yields a good approxima- 
tion to the standard Schmidt law with index 1.5. Less pre- 
dictably, the shock-induced rules also approximate Schmidt- 
type laws, albeit with somewhat shallower slopes. This 
outcome must depend on the dynamical structure of the 
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Figure CI. Global star formation in isolated disk simulations. 
Line types are identical to those in Fig. [2] in addition, the dot- 
dashed line is n = 1, m = 1, C* = 2.0. Note that here the top 
panel shows log M* . The simulations using shock-induced star 
formation produce brief initial bursts, but settle down promptly; 
no other features are seen in these plots. 

disks; for a hypothetical and rather extreme example, a 
bar-unstable disk would presumably yield different star for- 
mation profiles before and after the bar formed. A shock- 
induced model in which the average shock strength is in- 
dependent of R should yield a Schmidt-type law with an 
index of unity, regardless of the chosen value of m. Here, 
however, the shock-induced models with m — 1.0 appear to 
yield slightly steeper E, profiles than does the model with 
m = 0.5, hinting that the shocks are somewhat stronger near 
the center of this disk. 
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